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We examine sedimentation density profiles of star polymer solutions as an example of colloidal 
systems in sedimentation equilibrium which exhibit reentrant melting in their bulk phase diagram. 
Phase transitions between a fluid and a fluid with an intercalated solid are observed below a critical 
gravitational strength a*. Characteristics of the two fluid-solid interfaces in the density profiles 
occurring in Monte Carlo simulations for a < a* are in agreement with scaling laws put forth in the 
framework of a phenomenological theory. Furthermore we detect density oscillations at the fluid-gas 
interface at high altitudes for high gravitational fields, which are verified with density functional 
theory and should be observable in surface scattering experiments. 
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I. INTRODUCTION 



Colloidal particles in a suspension under gravitational 
influence show spatial inhomogeneities due to the symme- 
try breaking induced by the gravitational field. The prob- 
lem of sedimentation of particles in the presence of grav- 
ity has been of long scientific interest. The simplest ap- 
proximation is the one of non-interacting particles, valid 
in the limit of dilute solutions. This approach leads to 
an exponential sedimentation density profile, which was 
observed by Perrin for a calculation of Boltzmann's con- 
stant in 1910 Q|. Taking into account particle interac- 
tions at higher concentrations will yield corrections to the 
exponential density profile. For very small gravitational 
strength, a local-density-approximation (LDA) of den- 
sity functional theory (DFT) is justified ^ [|. In this 
case, there is a one-to-one correspondence between the 
sedimentation density profile and the isothermal equa- 
tion of state. This fact was exploited to extract the 
hard sphere equation of state experimentally by inves- 
tigating sterically stabilized colloids 0. Furthermore, 
within the LDA, a change in the height z corresponds to 
a local change of the chemical potential fj, of the bulk sys- 
tem. This implies that, in the limit of small gravity, the 
phase behavior becomes visible as a function of height z, 
a feature which was also been exploited to estimate the 
hard sphere freezing transition 0. Surprisingly, com- 
parison with Monte-Carlo (MC) simulations show that 
the LDA is even reliable for relatively strong inhomo- 
geneities or gravitational strengths This was further 
confirmed by comparing LDA against the exactly soluble 
hard rod model in one spatial dimension. While the LDA 
yields a monotonic decaying density profile p{z), a layer- 
ing shows up near the hard wall of the container bottom. 
Even crystallization can be induced by the bottom wall 
H . As shown recently || , details of this surface-induced 
crystallization may be significantly influenced by a pe- 
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riodic wall pattern. Indeed, pure colloidal crystals can 
be grown from sedimentation on a patterned substrate 
0, [§[ H| . In this case the gravitational field acts as an ex- 
ternal force enforcing and accelerating heterogenous nu- 
cleation and growth. Other fascinating phenomena in a 
gravitational field relevant for colloidal suspensions are 
phase transitions such as wetting Jlpj] , surface melting 
pd| , as well as dynamical effects as shock like fronts [[l2| , 
metastable phase formation jl3[l , long-range velocity cor- 
relations JlJ] , stratification |15f , and crystal growth jl6| . 

While the equilibrium sedimentation of hard sphere 
suspensions is well-understood g, g, g, |l7| jl§, charged 
suspensions are much more subtle as they reveal an ap- 
parent mass which is smaller than the bare mass at least 
for intermediate heights [g, [l9| 21. In this paper, 
we study a third kind of effective interaction between 
colloids, namely a very soft core as realized for star poly- 
mer solutions [p2| . The qualitative new feature of those 
solutions as compared to the traditional hard-sphere and 
Yukawa interactions is that their phase diagram exhibits 
a reentrant melting behavior for increasing density [ p3[ . 
In fact, our analysis holds for any system with a reentrant 
melting behavior but we will mainly focus explicitly on 
star polymers. Star polymers consist of / linear polymer 
arms attached to a central common core. The complete 
bulk phase diagram for star polymers in a good solvent 
was calculated in j2j| and exhibits several unusual solid 
lattices as well as reentrant melting. As will be discussed 
in detail in the following sections, due to the reentrant 
melting behavior, unusual density profiles, featuring in- 
teresting effects, arise and a wealth of scaling laws can 
be established. 

The paper is organized as follows: In Sec. II results of 
computer simulations of a system of star polymers, in- 
teracting by means of an ultrasoft pair potential pi| are 
presented. In Sec. we present a phenomenological 
theory giving account of the sedimentation profiles ob- 
served in the computer simulations. Scaling laws are put 
forth. Also in Sec. Ill, density functional theory in a sim- 
plified hybrid weighted density approximation (HWDA) 
is used to reproduce density oscillations at the fluid-gas 
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interface found in the simulation data, 
marks are contained in Sec. 



IV 



Concluding re- 



II. COMPUTER SIMULATION 

We performed canonical MC computer simulations 
keeping particle number TV, volume V , and temperature 
T constant. We used a simulation box with squared pe- 
riodic boundary conditions in x, y-direction and semi- 
infinite geometry in z-direction where the particles were 
confined only by the gravitational field for z > 0. The 
bottom wall at z = was hard and interacting with the 
star polymers by means of an effective star polymer-wall 
potential which is derived from the effective star polymer- 
hard sphere interaction in the limit of a sphere with zero 
curvature. The calculation was performed in |2J|. It is 
of the following form: 

pV sw (z) = A/ 3 / 2 

) z < 

&-m(f)-(^-l)(£i-i) *<f 

6(1 - erf(2Kz))/(l - erf(Ko-)) else. 

(1) 

With z we denote the distance from the center of one 
star polymer to the surface of the flat wall, a defines the 
so-called corona diameter of a star polymer, which is re- 
lated to its diameter of gyration cr g through a ~ 0.66cr g , 
see |2j|. The constants are A = 0.24, kg = 0.84, 

6 = 1/(1 + 2k 2 ct 2 ), & = ^iexp( K 2 ( 7 2 )(l-erf( K a)) 
and the inverse thermal energy (3 = 1/ksT. We empha- 
size that the range of the star-wall interaction is of the 
order of one or two corona diameters, so that the behav- 
ior of the sedimentation profiles for larger distances is not 
influenced. The star polymer pair potential is ultrasoft 
and is described by the following equation p4|: 



fiV ss {r) 




r < cr 



-V7/2 



a)) else, 



(2) 



with center-to-center distance r. Both interactions are 
purely entropic, hence they scale linearly with temper- 
ature. Previous work [^3| showed that a system of star 
polymers interacting by means of the potential pos- 
sess a very rich and interesting bulk phase diagram, see 
Fig. 0, exhibiting reentrant melting and reentrant freez- 
ing transitions for arm numbers f c < / < 54, with the 
critical arm number f c = 34. As we will discuss in more 
detail below, it is the reentrant melting that makes this 
type of system appropriate for the analysis presented in 
this article. The suspending liquid is assumed to be in- 
compressible. Furthermore we treat the solvent to be 
continuous, neglecting possible effects of the discreteness 
of the solvent particles. Given the size of the colloidal 
particles under observation, the star polymers, this is 



a reasonable assumption. In the simulation, the initial 
configuration of the system was chosen to be a body- 
centered cubic (bcc) solid to facilitate equilibration. Its 
lattice constant a was determined from a bulk system 
with a packing fraction 77 = ^pcr 3 — 0.5 lying in the bcc- 
regime in the bulk phase diagram, see Fig. [j]. The lateral 
box-dimensions were chosen to be multiples of the lattice 
constant a. 
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FIG. 1: Bulk phase diagram of star polymers interact- 
ing with potential (^|), calculated in p3[ . Arm num- 
ber / is plotted versus packing fraction r/. The squares 
indicate the phase boundaries; solid lines are guide to 
the eye. The black cross denotes the point with crit- 
ical arm number f c ~ 34 and corresponding density 
?y c ~ 0.43. The system is always fluid for arm numbers 
smaller than the critical arm number f c and shows reen- 
trant melting behavior for arm numbers f c < / < 54. 
The arrow indicates a path through the phase diagram 
that is equivalent with a change in the altitude z within 
the LDA for small a. The four observed solid phases 
are body-centered cubic (bcc), face-centered cubic (foe), 
body-centered orthogonal (bco), and diamond (diam.). 



The total number of particles was then fixed by prescrib- 
ing a certain value of the thermodynamic variable r, giv- 
ing the number density per unit surface. The density 
profile p(z) is normalized as 



p(z)dz. 



(3) 



to 2 is the number of particles piled up over the area a 2 
of the bottom wall. Typical system sizes were N = 2000 
particles and the Monte Carlo runs were extended over 
JVmc ~ 500 000 cycles, each cycle comprising one trial 
move for each of the N particles. Besides the aforemen- 
tioned thermodynamic variable r, two further parameters 
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characterize the state of the system: First, the arm num- 
ber / of the star polymers, being the number of polymer 
chains grafted on the central core. Second, the dimen- 
sionless gravitational strength (or Peclet number) 



mga 



(4) 



which describes the ratio of the potential energy gain to 
the thermal energy ksT for a particle of mass m, dis- 
placed by a in height in an external held with acceler- 
ation g. The three parameters /, r and a were varied 
over a broad range of values. The particles were moved 
by employing the standard Metropolis algorithm. 



In Fig. U we show results for different gravitational 
strengths a, while to 1 — 48.87 and / — 39 are fixed. 
The / = 39 star polymer system displays reentrant melt- 
ing in the fluid — > bcc — > fluid sequence, as seen along 
the arrow in Fig. |[ The gravitational field forces the 
local density cr 3 p(z) to take values that scan the range 
from cr 3 p(z) = up to high values, cr 3 p(z) = 3. Thus, 
the local density 'crosses through' the range of the phase 
diagram where the system displays a bulk bcc phase. It 
it intuitively expected that the system will then feature a 
solid regime (for intermediate densities) intercalated be- 
tween two fluid regimes, at low and high densities. We 
have found that this is indeed what happens but provided 
that the gravitational strength does not exceed a critical 
value a* , as we discuss below. 

Let us start from the case where no solid phase ap- 
pears. For a > a* [Figs. ||(a) and (b)], we obtain den- 
sity profiles p{z) that show three distinct features: First, 
there is layering on the wall due to packing effects, typi- 
cally extending over several layers. As z increases a fluid 
regime with density decaying as a linear function of al- 
titude z can be distinguished. At some height (z ~ 25a 
in (a)) the density rapidly decays to zero. At this strong 
inhomogeneity oscillations in density with wavelength a 
can be distinguished in the sedimentation profile, which 
is smooth elsewhere in the linear regime. The linear de- 
pendence of the density profile on z, can be understood in 
terms of a local den sity fu nctional mean-field theory, as 
will be shown in Sec. HI B; the corresponding results from 
this theory are shown in Figs. ||(a) and (b) with dotted 
lines. The density oscillations observed in the simulations 
were reproducible in the framework of density functional 
theory using a simplified form of the H WDA, as will be 
discussed in further detail in Sec. 
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By lowering the gravitational strength a further, a crit- 
ical strength a* in the range 16.0 < a* < 17.0 is discov- 
ered. Below a* the density profiles qualitatively change 
and exhibit a new feature. Strong density oscillations 
appear, a clear indication for a crystalline phase. They 
extend over 10 to 20 star diameters, equivalent to several 
crystalline layers. The length of the crystal grows, as a 
decreases. A typical simulation snapshot is shown in Fig. 
|3| next to the corresponding equilibrium density profile. 
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FIG. 2: Sedimentation profiles of star polymers for an 
arm number / = 39 and a density to 2 — 48.87. The 
gravitational strength a is decreased from (a) to (f) 
with (a) a = 30.0, (b) a = 17.0, (c) a = 16.0, (d) 
a = 8.0, (e) a = 6.0 and (f) a = 4.0. In plots (c) - 
(f) the order parameter $4 is also shown (dashed line) 
using the same y-scale as the profiles. In (a) and (b) a 
straight line whose equation is derived within the LDA 
[see Eq. (12h] is superimposed on the plots (dotted line). 



Here, the well-ordered crystal phase in the middle of the 
simulation box (20c < z < 30cr) is clearly visible. 

As an additional check for crystalline order, we calcu- 
late the local order parameter, ^4, that checks for four- 
fold symmetry in two dimensions around a given particle. 
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It is defined by 



tf 4 (z) = 



1 Nl 

W, zJ zJ 



i=i <fe> 



(5) 



where the fc-sum includes the four nearest neighbors of 
the given particle and the j-sum extends over N\ particles 
in the corresponding layer. A layer is defined by a slab of 
thickness 5 ~ 0.2 a, centered around the given particle at 
elevation z, which is motivated by the 'Lindcmann melt- 
ing rule', assuming a maximum particle displacement of 
approximately 10% around the equilibrium position in a 
possible crystal regime. The angular brackets indicate 
a canonical ensemble average. <j)jk is the polar angle of 
the interparticle distance vector with respect to a fixed 
reference frame. For ideal fourfold symmetry, i.e., for a 
particle contained in a bec-solid layer, $4 = 1. Due to 
thermal motion, small defects of the perfect crystalline 
symmetry arise and usually values of \& 4 > 0.8 || are 
taken to be conclusive evidence for a crystalline phase 
with fourfold-in-layer-symmetry. As can be seen in Figs. 
||(c)-(f) our simulation data do indeed show values up to 
^4 1=3 0.95 in the region of the density profile p(z) which 
we already identified to be solid due to the pronounced 
density oscillations. 

Comparing the interval of the packing fraction in which 
crystallization occurs to the bulk phase diagram in Fig. [j], 
we may thus conclude that the intercalated solid regime 
is a manifestation of the reentrant melting in the bulk 
phase diagram, mapped onto the z-axis in a system un- 
der gravitational influence. The absence of freezing for 
strong gravitational fields (a > a*) can now be at least 
qualitatively understood: for high values of a, the den- 
sity profiles grow too fast as z approaches the wall, so 
that the mapping onto the z-axis results into a domain 
which is too narrow to sustain cr ystall ine order. In fact, 
as we will show in detail in Sec. IIIB, a minimal, non- 



vanishing thickness of the crystalline layer is necessary 
so that the latter can be stably 'nested' between the two 
fluid phases. 




FIG. 3: Snapshot from MC simulation (right) shown with 
the corresponding equilibrium density profile (left). The 
star polymers are rendered as spheres with diameter a. 
The parameters are: / = 39, a = 5.0, to 1 — 42.1. 
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III. THEORY 

A. Density Functional Theory in Local Density 
Approximation (LDA) 

In order to predict scaling relations characterizing crys- 
tallization in sedimentation profiles of star polymer solu- 
tion we apply density functional theory within the frame- 
work of the local density approximation (LDA) . The lat- 
ter is a reliable theoretical tool in cases where the density 
profile of the system varies slowly with z, so that it can 
be considered as staying essentially constant at length 
scales set by the microscopic natural length of the sys- 
tem (cr in this case). As can be seen in Figs. ||(a) and 
(b), this is indeed the case if we discard the strong os- 
cillations close to the wall (the layering effect). As the 
range of these oscillations is much shorter than the range 
of the density profile itself, the bulk of the free energy of 
the system resides in the smooth "ramp-like" part of the 
density profile and the use of the LDA is justified. Ac- 
cordingly, we will omit the star-wall potential from our 
considerations in this subsection and consider only the 
effects of the external gravitational field $ C xt( z ) = mgz. 

We work in the grand canonical ensemble and intro- 
duce the chemical potential p and a variational grand po- 
tential per unit area, E(T, p; |jo(z)]) which is a functional 
of the density profile. Introducing the ideal and excess 
per unit area contributions to the intrinsic Helmholtz free 
energy of the system, F i( i[p(z)] and F cx [p(z)] respectively, 
we find that in the LDA, the expression for E(T, p; [p(z)]) 
reads as: 

E(T, /*,[/>(*)]) = F id [p(z)] + F ex [p(z)} 

+ / dz$ cxt (z)p(z) -pi dzp(z) 



k B T / Azp(z) [ln(p(z)A 3 ) - l] 
Jo 



&z [f(p(z)) + {mgz - p)p{z)} , 



(6) 



where A = y/hP^hrmk^T is the thermal de Broglie wave- 
length and f(p{z)) is the Helmholtz free energy per unit 
volume of the bulk fluid. The minimization of E with 
respect to p{z) yields the equilibrium profile po(z); the 
value of the functional at equilibrium, E(T,/i, [po^)]) is 
then the Gibbs free energy per unit area, E(T, p) of the 
system. Setting St(T, p,[p(z)])/Sp(z)\ pn{z) = in Eq. 
(0), leads to: 



k B T\n(p (z)a 3 ) + f'(p (z)) = p! - mgz, 



(7) 



where f'(x) denotes the derivative of f(x) and p! = p — 
31n(A/<r) is a shifted chemical potential. 

Due to the ultrasoft character of the logarithmic- 
Yukawa star-star interaction V ss (r), the star polymer sys- 
tem belongs to the class of mean-field fluids | 26 , ^7], |2^, 



p9| , pjQ |, for which the excess free energy density is a 
quadratic function of p, namely: 

f(p) = Y Air J dr r Vss ( r > =: — 2 — ' ( ' 

with the Fourier transform V ss (k) of the pair potential. 
This property is valid for high density fluids provided 
their pair potential V(r) is only slowly diverging at the 
origin and decays fast enough to zero as r — > oo, so 
that it is integrable. For the more restrictive case of 
a nondiverging potential at r = 0, the stronger condition 
c(r) = -/3V{r) holds approximately @ H, ll> HI' with 
c(r) denoting the direct correlation function of the fluid 
pl| . This gives again rise to Eq. (| ) above through the 
compressibility equation of state (31 ] . 

Using the dimensionless variables x = z/cr, p{x) = 
p(z)a 3 , B = (3V ss (0)/a 3 and p = f3p' and introducing 
Eq. (||) into Eq. M), we obtain the equilibrium profile 
through the equation: 



\n(p (x)) + Bp (x) — p — ax. 



(9) 



For star functionality / = 39 we obtain B = 250 and, for 
/ = 32, B = 204 [see Eq. (§)]. Hence, the second term in 
the Ihs of Eq. (JsJ) above dominates over the logarithmic 
term for densities p(x) > 0.10. As almost the entire 
simulation density profile fulfills this condition, we finally 
omit the logarithmic term from Eq. (|^) above and obtain 
thereby a linear density profile: 



Po(x) 



for x < 0, 
— g — for < x < p/a, 
for p/a < x. 





ft— ax 



(10) 



The chemical potential p is now determined through the 



normalization condition J Q 
ing: 



&xpq(x) — to 2 = r, yield- 



= V2aBf, 



(11) 



and from Eq. @ the final expression for the density 
profile: 



Po{x) 







for x < 0, 
\x for < x < y^f? 
for ,[W < x. 



(12) 



The prediction (12) is compared against the MC simu- 
lation results in Figs. ||(a) and (b); theory and simulation 
are in excellent agreement. This linear dependence of the 
density profile on z is the first scaling prediction we make 
for such systems. Moreover, by introducing Eq. (10) into 
Eq. (|^) , and once more ignoring the logarithmic term, we 
find that the Gibbs free energy per unit area S(T, p) is 
a power-law of the chemical potential, namely: 



/3a 2 E(T,/2) = 



6aB 



(13) 
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Accordingly, the Helmholtz free energy per unit area, 
(3a 2 F(T, f) = /3cr 2 S(T, p) + pf obeys the scaling law: 



(3a 2 F(T,?) = -V2aB? 3 / 2 . 



(14) 



The thermodynamic relation p, = d(Pa 2 F)/dr returns 
Eq. @. ' 

We now examine whether the density oscillations oc- 
curring at high z-elevations, which are clearly visible in 
Figs. |2|(a), (b), can be obtained in the framework of the 
full LDA, with the logarithmic term included, Eq. 
Though the latter is an implicit equation for p{x) , we do 
not need to solve it in order to answer the question at 
hand. The key observation is that Eq. @ delivers an 
explicit functional form for the inverse function: 



x(p) = — a 1 (In p + Bp — (i). 



(15) 



If the LDA profile displayed oscillations, then p(x) would 
go through various maxima and minima and there should 
be several points x m where the derivative p'(x m ) would 
vanish, with the implication that the derivative of the in- 
verse function, x'(p m ), would diverge at the correspond- 
ing density values p m . From Eq. (Ilq) above, we obtain 
x'(p) = -a^ 1 (p^ 1 + B) < for all < p < oo. The only 
divergence of x'(p) occurs for the trivial limit p — > and 
corresponds to the exponential decay p{x) cx e~~ ax , valid 
for high elevations. The LDA is incapable to reproduce 
this effect, a feat that, in fact, could have been antici- 
pated: these oscillations occur at length scales a, whereas 
the LDA is applicable when the spatial inhomogeneity of 
the profile has a ch aract eristic length much larger than 
the latter. In Sec. Ill C , we resort to a more powerful 
density functional approximation in order to reproduce 
this feature of the density profile. 



B. Phenomenological Landau Theory 

If one focuses close to the reentrant melting transition 
point and studies the case for small gravity (a -C 1), a 
phenomenological Landau-like approach can be adopted 
to explore further scaling predictions for the crystalliza- 
tion transition. We study the situation sketched in Fig. 
^ of a crystalline sheet of width I intervening between two 
fluid parts of the sedimentation profile. Let us define the 
excess grand canonical free energy per unit area, S ox (Z), 
in such a situation with respect to a situation where no 
crystallization takes place. This quantity is given by: 



£«(/) = E(0-S(i = 0), 



(16) 



with being the grand canonical free energy per unit 
area when a solid of thickness I is present. Evidently, 
£(Z = 0) is the quantity given by Eq. ([l3[) above. 

The excess grand canonical free energy S ex per unit 
area comprises of three parts: 

1. The equilibrium surface tensions 71 and 72: These 
describe the additional free energy cost in creating 
the two solid- fluid interfaces at z = Z\ and z = z%. 




FIG. 4: Sketch of the situation in which an in- 
tervening solid of width I is nested between two 
fluids under the influence of a gravitational field. 



2. A thermodynamic contribution which essentially 
depends on the arm number. If / > f c this contri- 
bution favors a solid sheet. 

3. A free energy penalty due to an elastic distortion 
of the solid in the external field. 



Hence: 



S cx = 71+72 



-'clast ■ 



(17) 



Let us discuss the different contributions in more detail: 
The surface tension will mainly control the relative ori- 
entation of the solid with respect to the z-direction. One 
expects that a close-packed surface of the bec-solid (i.e. 
a (100) orientation) will have smallest surface tension 
and will hence be the realized orientation. In fact, this 
is what we found in our simulation data presented in 
Figures ^ and ||. For hard-sphere fcc-solids the interfa- 
cial fluid-solid free energy has been calculated recently 
in equilibrium by computer simulation p^| . Its order of 
magnitude is 



, (* = 1,2), 



(18) 



where a is a microscopic length scale. 

The thermodynamic contribution could be calculated 
within the LDA with f(p, T) taken from liquid state the- 
ories for the fluid and solid cell theory for the crystal. 
Here, we will simply focus on a Landau-type theory close 
to the reentrant melting point characterized by a criti- 
cal star number density p c = -^ji]c and the critical arm 
number / c , see Fig. |l|. Performing a Landau expansion 
and dropping the temperature dependence one gets 



Hp) = fs(Pc) + M&-f) 

+f' S (Po)(p- Pc) 

+ lfs(Pc)(p-Pc) 2 - 



(19) 



Here, f s {p) is the free energy per unit volume of the solid 
phase and A s is a constant governing the first leading 
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term in an expansion around / = / c . Likewise in the 
fluid phase one has 



Hp) = ff(p c ) + Mf c -f) 

+/f(pc)(p-pc) 



(20) 



with / s (p c ) = f f {p c ), f s { Pc ) = / f '(p c ), but / 8 "(p c ) > 
ff(Pc) m general. Performing the inversion of /(p) in 
order to get the density profile leads to a piecewise linear 
profile for the averaged density with two density jumps 
at z = Z\ and z — Z2, as determined by the Maxwell 
construction, see Fig. |J. 



+ p c for zi < z < z 2 

/ f "(Pc) elSe - 



(21) 



Consequently, by inserting this into the free energy func- 
tion one gets 



Std = -a(f - f c )l + 



1 



f 



/f"(Pc) / s "(Pc) 



m 2 g 2 l 3 



12 



(22) 



Note that a > in order to stabilize the solid for / > / c . 

Third, the elastic part can be calculated by contin- 
uum elastic distortion theory of the solid. For a different 
situation of a solid in an external field this has been for- 
mulated by Gittes and Schick |33| . Following these ideas, 
we assume a z-independent lateral strain ey but consider 
a z-dependent vertical strain e±. By symmetry, ey has to 
be zero for the crystal being stable at z = zq, i.e., for the 
crystal at the reentrant melting point. Elasticity theory 
predicts for S i ast 



1. The realized crystalline thickness I as obtained by 
minimizing E ox with respect to I for fixed a and / 
scales as 



I 



(27) 



2. The phase transition to a crystalline sheet is first 
order. It happens beyond a critical a-dependent 
arm number / cr jt where 



/crit - fc OC 7 ' OL 



2/3_2/3 



(28) 



with a scaling exponent of |. 

3. The width / corresponding to the transition scales 
as 



i oc 7 1/3 «- 2/a oc-^— 

J Jc 



(29) 



The analysis presented here is general, the scaling pre- 
dictions derived are valid for any reentrant melting be- 
havior in equilibrium (e.g., laser-induced freezing |fi4| , 
or polydisperse systems Furthermore all these re- 

lations can, in principle, be checked by simulation. Rela- 
tions ( p8| ) and (29), however, require high computational 
efforts. In order to check on scaling relation (|27|), we 
measured the crystal length I in MC simulations varying 
a or /, while keeping the density r fixed. The crystal 
length is determined by the range Az, where the order 
parameter ^(z) has values larger than 0.8. The results 
are plotted in Fig. |5|, showing excellent agreement with 
the scaling predictions. 



(23) 



Weighted Density Approximation of the 
Density Functional 



where C > is related to the elastic constants of the 
solid. As 



we obtain 



e_L oc p — p c oc mgz, 



Seiast = C'l 3 m 2 g 2 , 



(24) 



(25) 



with another constant C . Eq. (|25|) has a similar form as 
the second term of Eq. (^2|). 

In summary, the total grand canonical excess free en- 
ergy is 

£«(Q = -a(f - f c )l + ba 2 l 3 + 7l + 72 , (26) 

where a, b > 0. Moreover we will use 7 = 71 + 72 from 
now on. A first order phase transition takes place if the 
minimum of S cx with respect to I yields S ex = 0. 

This determines the following resulting scaling rela- 
tions: 



In order to verify the density oscillations close to the 
fluid-gas interface of the sedimentation profiles for large 
values of a, we apply a simplified form of the HWDA (hy- 
brid weighted density approximation). The full HWDA 
was constructed by Leidl and Wagner in |37). Given an 
external field $ C xt{z) the free energy is a unique func- 
tional of the density profile p(z). Thus, the excess free 
energy per unit surface in the HWDA framework is given 
by 



T 



Ap\ = f 

Jo 



p(z)f (p(z))dz, 



(30) 



where fo(p) denotes the excess free energy per particle 
of a homogeneous liquid of density p. The weighted den- 
sity p(z) follows from a convolution with the weighting 
function u)(r; p) 



p(z)= / p(z>(|r-r'| ;/ 3)dr' 



(31) 
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FIG. 5: Verificatio n of t he scaling behavior theoretically 
predicted in Sec. [[II B in Eq. (B7|), calculated by MC 
simulations. In (a) the crystal length I is plotted versus 
V f — fc, keeping a — 6.0 fixed. In (b), I is plotted ver- 
sus 1/a for a fixed arm number / = 39. The dashed 
lines are linear fits to the simulation results (circles). 



with a global density p. The weighting function w[r;p) 
is fixed by a simple quadratic equation in Fourier space 



2/o(Po)w(fc;p ) 



Po/o (po)u 2 (k;po) 



(32) 



The primes denote differentiations with respect to the 
density p and c^(k;po) is the Fourier transform of the 
direct correlation function of the homogeneous fluid. A 
unique solution of u>(k; po) is determined by the normal- 
ization Q(k = 0; po) = lj also ensuring the compressibil- 
ity rule to be satisfied. We have solved the homogeneous 
problem with Ornstein-Zernike fluid integral equations 
using the Rogers- Young closure Q. Resulting corre- 
lation functions and structure factors are in very good 
agreement with MC simulations of the bulk system [^6| . 
In difference to the complete HWDA, where the global 
density p is chosen to be a functional of p[z), we keep 
p fixed. This simplification is sufficient to verify the ob- 



served oscillations, accompanied by the advantage that 
the numerical effort is enormously reduced. Best agree- 
ment with simulation results could be achieved when 
choosing the global density p to be of the order of the 
averaged density near the bottom wall z = 0. The tails 
of the density profiles are nearly unaffected by the choice 
of p. A similar approach is used in the SWDA (3{| for 
an inhomogeneous fluid in contact with a bulk fluid of 
density pb] there p was chosen to be pb- 

Applying the usual Euler-Lagrange minimization for 
the Helmholtz free energy per unit area T\p\ with chem- 
ical potential p 



5T 



Sp(z) 



= M - ^ext(^)) 



(33) 



and using <& cxt (z) = 
density profile p(z) 



Vsw(z) + az/ (3a, we obtain for the 



Zexp{c<-V(z; [p]) - az/a - (3V sw (z)} z>0 



(34) 







else. 



The fugacity £ is determined by the normalization con- 
dition £ = rj J °° dz exp{c^' (z; [p]) — az/a — (3V sw (z)}. 
c^(z; [p]) is the one-particle correlation function: 



-i c (D 



Sp(z) 



(35) 



+j dz'p(z')ti(p(z')M\z-z'\;p). 

Equation (Q) was solved for the profile by standard it- 
erative techniques, see, e.g., Ref. Q. The results for an 
arm number / = 32, ra 2 = 21.8 and three different val- 
ues of a are shown in Fig. || together with MC simulation 
data. The global density p for all three profiles is fixed at 
pa 3 = 1.8. The DFT results are in very good agreement 
with the simulation profiles. In particular, the interface 
oscillations with wavelength a also occur in the DFT. For 
a = 10.0, the lowest a that is shown, the profile is nearly 
indistinguishable from MC data far from the wall, while 
for increased a, a — 20.0 and a = 30.0 the interface 
oscillations are underestimated in our DFT approach. 

These oscillations are not a specific feature of the star 
polymers; we have also performed MC simulations using 
a repulsive Yukawa interaction of the form 



V(r) oc exp(— Kr)/r, 



(36) 



corresponding to the part of the star polymer pair po- 
tential (Q) valid for distances r > a. Here exactly the 
same behavior could be found at the tails of the density 
profiles. For hard sphere systems, on the other hand, 
no such density oscillations are present. This suprising 
fact might be attributed to the long ranged tail in the 
interaction potential. 
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FIG. 6: Density profiles p(z)a 3 for three different values 
of a = 10.0, 20.0, 30.0 and fixed arm number / = 32 
and fixed density rer 2 = 21.8 calculated with DFT 
(solid curves) compared to MC simulation results (dashed 
curves). The slope of the curves increases with in- 
creasing a. The inset shows the bare DFT results 
for a better identification of the interface oscillations. 

IV. CONCLUSIONS 

Concluding, we have presented results for systems ex- 
hibiting reentrant melting in the bulk phase diagram, un- 



der gravitational influence. It was shown that a phase 
transition occurs when the gravitational strength a is 
varied: Below a critical a*(f,r), intercalated crystal- 
lization occurs in the sedimentation profiles of the ob- 
served star-polymer solutions, whereas for a > a* (/, t) 
we find monotonic sedimentation profiles p{z). In MC 
computer simulations scaling relations for the crystalliza- 
tion, predicted in the framework of a phenomenological 
theory, valid for all systems exhibiting reentrant melting 
in the bulk phase diagram, could be verified. Using den- 
sity functional theory, density oscillations at the fluid-gas 
boundary, observed in the MC simulations, could be re- 
produced. 



In principle, our results can be verified in surface- 
sensitive scattering experiments or real-space imag- 
ing methods for colloidal suspensions. Unlike non- 
monotonicities on the liquid side of the gas- liquid equilib- 
rium interface (see e.g. Q), the density oscillations on 
top of the sedimentation profile are not affected by cap- 
illary wave fluctuations and may thus be verified in real 
samples. The intervening solid sheet should be signalled 
by a Bragg- like peak in surface reflection measurements. 
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